home *** CD-ROM | disk | FTP | other *** search
/ Meeting Pearls 1 / Meeting Pearls Vol 1 (1994).iso / installed_progs / gfx / lise2.1 / lise / src / augcal.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-31  |  8.0 KB  |  344 lines
  1. /*
  2.    calibrate energy scale of auger spectrum
  3. */
  4.  
  5. #define ECALTMP "ecal.auger"
  6.  
  7. #include <stdio.h>
  8. #ifdef AMIGA
  9. #include <gfxamiga.h>
  10. #else
  11. #include <gfx.h>
  12. #endif
  13.  
  14. float x,y,*spc,*err,*tim;
  15. int   fdlockin, delay;
  16. float stepv, startv, stopv;
  17. float ecal;
  18. int spcmax;
  19. char *cname[110];
  20. char *lname[110];
  21. float *Eauger, *Ebind, *sens, *Ratom;
  22.  
  23.  
  24. int find_peak(Epos)
  25. float Epos;
  26. {
  27. float    pmin, pmax, x, y, mean, diff;
  28. int      k, ka, kmin, kmax, isl, islm, isr, isrm;
  29.  
  30. /* search minimum */
  31.    ka = (int) (Epos / _tica);
  32.    isl = (int) ((5.0 / _tica) / 2.0);
  33.    islm = (int) (5.0 / _tica);
  34.    isr = (int) ((5.0 / _tica) / 2.0);
  35.    isrm = (int) (5.0 / _tica);
  36.  
  37.    while((isl < islm) && (isr < isrm)) {
  38.       if(((ka - isl) < 1) || ((ka + isr) > spcmax)) return(0);
  39.       pmin = spc[ka];
  40.       kmin = ka; mean = 0.0;
  41.       for(k = (ka - isl); k <= (ka + isr); k++) {
  42.          mean = mean + spc[k];
  43.          if(spc[k] < pmin) {
  44.             pmin = spc[k];
  45.             kmin = k;
  46.          }
  47.       }
  48.       mean = mean / ((float) (isr + isl + 1));
  49.       if(spc[ka - isl - 1] <= pmin) {
  50.          isl = isl + 1;
  51.          continue;
  52.       }
  53.       if(spc[ka + isr + 1] <= pmin) {
  54.          isr = isr + 1;
  55.          continue;
  56.       }
  57.       diff = ((mean - pmin) / mean); if(diff < 0.0) diff = -1.0 * diff;
  58.       if(diff < 0.01) {
  59.          isr = isr + 1;
  60.          isl = isl + 1;
  61.          continue;
  62.       }
  63.       break;
  64.    }
  65.  
  66.    diff = ((mean - pmin) / mean); if(diff < 0.0) diff = -1.0 * diff;
  67.    if(diff < 0.01) return(0);
  68.  
  69. /* search maximum */
  70.    isl = (int) (8.0 / _tica);
  71.    islm = (int) (12.0 / _tica);
  72.  
  73.    kmax = kmin;
  74.    pmax = pmin;
  75.    while(isl < islm) {
  76.       if((kmax - isl) < 1) return(0);
  77.       for(k = (kmax - isl); k <= kmax; k++) {
  78.          if(spc[k] > pmax) {
  79.             pmax = spc[k];
  80.             kmax = k;
  81.          }
  82.       }
  83.       if(spc[kmin - isl - 1] >= pmax) {
  84.          isl = isl + 1;
  85.          continue;
  86.       }
  87.       break;
  88.    }
  89.    x = _tica * ((float) kmin);
  90.    y = _tica * ((float) kmax);
  91.  
  92. /* Difference between maximum and minimum should not exceed 20 eV*/
  93.    diff = x - y ; if(diff < 0.0) diff = -1.0 * diff;
  94.    if(diff > 20.0) return(0);
  95.  
  96. /* Difference between measured peak and value of table should not exceed 20% */
  97.    diff = x - Epos; if(diff < 0.0) diff = -1.0 * diff;
  98.    if((diff / Epos) > 0.2) return(0);
  99.  
  100.    if(kmax == kmin) return(0);
  101.  
  102.    return(kmin);
  103. }
  104.  
  105.  
  106. float find_energy(z)
  107. char *z;
  108. {
  109. int i,n,m;
  110.  
  111.    z[0] = z[0] & 223;               /* = toupper ! */
  112.    if(z[1] > 0) z[1] = z[1] | 32;   /* = tolower ! */
  113.    for(n = 0; n < 110; n++) {
  114.       if(strcmp(cname[n],z) == 0) break;
  115.    }
  116.    if(n >= 109) {
  117.       printf("sorry, element >%s< not found\n");
  118.       exit(0);
  119.    }
  120.    return(Eauger[n]);
  121. }
  122.  
  123.  
  124. help()
  125. {
  126.   printf("augcal spectrumfile [-chan n] [-cur] -energy n.m\n");
  127.   printf("       change energy calibration of specified file\n");
  128.   printf("  -chan n           specify channel\n");
  129.   printf("  -cur              use cursor to select the channel\n");
  130.   printf("  -energy n.m       set energy of peak\n");
  131.   printf("  -line Xx          use auger energy of element Xx as energy\n");
  132.   printf("  -auto Xx          perform auto calibration for element Xx\n");
  133.   printf("  -table filename   uses a custom table for elements\n");
  134.   printf("  -v                verbouse: print out calibration factor\n");
  135. }
  136.  
  137. float frline(stream)
  138. FILE *stream;
  139. {
  140. int i,n;
  141. char c,s[80];
  142.  
  143.    fgets(s,80,stream);
  144.    return(atosf(s));
  145. }
  146.  
  147. main(argc,argv)
  148. int argc;
  149. char *argv[];
  150. {
  151. int n,m,i;
  152. char *z,*comment, *elements;
  153. int channel;
  154. int flag_v;
  155. float energy;
  156. float real_t,real_y;
  157. FILE *fp;
  158.  
  159.    flag_v = 0;
  160.    energy = 1.0;
  161.    channel = 1;
  162.  
  163.    if(argc < 4) {help(); exit(0);}
  164.  
  165.    fdlockin = auxopen("auger");
  166.    delay = atoi(auxparams[4]);
  167.    stepv = atosf(auxparams[5]);
  168.    startv = atosf(auxparams[6]);
  169.    stopv = atosf(auxparams[7]);
  170.    ecal = atosf(auxparams[8]);
  171.    close(fdlockin);
  172.  
  173.    Eauger = (float *) calloc(110,sizeof(float));
  174.    Ebind = (float *) calloc(110,sizeof(float));
  175.    sens = (float *) calloc(110,sizeof(float));
  176.    Ratom = (float *) calloc(110,sizeof(float));
  177.    z = (char *) malloc(80);
  178.    comment = (char *) malloc(80);
  179.    elements = (char *) malloc(80); strcpy(elements,"elements");
  180.    spc= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  181.    err= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  182.    tim= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  183.    if(tim==NULL) {
  184.       printf("sorry, not enough memory\n");
  185.       exit(-1);
  186.    }
  187.  
  188.    for(n = 0; n < 110; n++) {
  189.       cname[n] = (char *) malloc(4); cname[n][0] = 0;
  190.       lname[n] = (char *) malloc(16); lname[n][0] = 0;
  191.       sens[n] = 0.0;
  192.       Eauger[n] = 0.0;
  193.       Ebind[n] = 0.0;
  194.       Ratom[n] = 0.0;
  195.    }
  196.  
  197.    spcmax=readspec(argv[1],spc,err,tim,comment);
  198.  
  199.    if(checkopt(argc,argv,"-v",z)) flag_v = 1;
  200.    if(checkopt(argc,argv,"-chan",z)) channel = atoi(z);
  201.    if(checkopt(argc,argv,"-cur",z)) {
  202.       getcur(&channel,&real_t,&real_y);
  203.       if(flag_v == 1) printf("channel = %d\n",channel);
  204.    }
  205.    if(checkopt(argc,argv,"-energy",z)) energy = atosf(z);
  206.    if(checkopt(argc,argv,"-table",z)) strcpy(elements,z);
  207.    if(checkopt(argc,argv,"-line",z)) {
  208.       read_elements(elements);
  209.       energy = find_energy(z);
  210.    }
  211.    if(checkopt(argc,argv,"-auto",z)) {
  212.       read_elements(elements);
  213.       energy = find_energy(z);
  214.       channel = find_peak(energy);
  215.       if(channel == 0) {
  216.          printf("could not find peak !\nNo calibration performed !\n");
  217.          exit(0);
  218.       }
  219.    }
  220.  
  221.    _tica = energy / channel;
  222.    for(n = 0; n <= spcmax; n++) tim[n] = _tica * ((float) n);
  223.    writespec(argv[1],spc,err,spcmax,2,comment);
  224.    strcpy(z,argv[1]); strcat(z,".tim");
  225.    writespec(z,tim,err,spcmax,2,comment);
  226.  
  227.    ecal = _tica / stepv;
  228.    if(flag_v == 1) printf("calibration is %f eV/V\n",ecal);
  229.  
  230.    fp = fopen(ECALTMP,"w");
  231.    if(fp != NULL) {
  232.       fprintf(fp,"%f\n",ecal);
  233.    }
  234.    free(err); free(spc); free(tim);
  235.    free(elements); free(comment); free(z);
  236.    for(n = 0; n < 110; n++) {
  237.       free(lname[n]);
  238.       free(cname[n]);
  239.    }
  240.    free(Ratom); free(sens); free(Ebind); free(Eauger);
  241.    exit(0);
  242. }
  243.  
  244. getcur(channel,time,value)
  245. int *channel;
  246. float *time, *value;
  247. {
  248. float tmin,
  249.       tmax,
  250.       ymin,
  251.       ymax,
  252.       real_t,
  253.       real_y,
  254.       yfak,
  255.       tfak,
  256.       xoffset,
  257.       yoffset;
  258. int   leftmargin,
  259.       bottommargin,
  260.       rightmargin,
  261.       topmargin,
  262.       _win_flg;
  263. int x,xx,y,yy;
  264. FILE *fp;
  265.  
  266.       tekopen();
  267.       getxy(&x,&y);
  268.       close(_tek4014);
  269.       /* calculating channel and time equivalents */
  270.       fp=fopen(ASHBDRY,"r");
  271.       leftmargin = frline(fp);
  272.       bottommargin = frline(fp);
  273.       rightmargin = frline(fp);
  274.       topmargin = frline(fp);
  275.       _win_flg = frline(fp);
  276.       tmin = frline(fp);
  277.       tmax = frline(fp);
  278.       ymin = frline(fp);
  279.       ymax = frline(fp);
  280.       _tica = frline(fp);
  281.       fclose(fp);
  282.       yfak=(topmargin-bottommargin)/(ymax-ymin);
  283.       tfak=(rightmargin-leftmargin)/(tmax-tmin);
  284.       xoffset= ((- tmin) * tfak) + leftmargin;
  285.       yoffset= ((- ymin) * yfak) + bottommargin;
  286.  
  287.       real_t= (x - xoffset)/tfak;
  288.       real_y= (y - yoffset)/yfak;
  289.       xx = real_t / _tica;
  290.       *channel = xx;
  291.       *time = real_t;
  292.       *value = real_y;
  293. }
  294.  
  295. read_elements(elements)
  296. char *elements;
  297. {
  298. int i,n,m;
  299. char z[80];
  300. FILE *fp;
  301.  
  302.    fp = fopen(elements,"r");
  303.    if(fp == NULL) {
  304.       fprintf(stderr,"file >%s< not found !\n",elements);
  305.       exit(0);
  306.    }
  307.    while(!feof(fp)) {
  308.       fgets2(z,80,fp);
  309.       if(feof(fp)) break;
  310.       if(z[0] == ';') continue;
  311.       if(z[0] < 33) continue;
  312.       n = atoi(z);
  313.       if((n < 1) || (n > 108)) {
  314.          fprintf(stderr,"error in >%s< unexpected >%s<\n",elements,z);
  315.          exit(0);
  316.       }
  317.       fgets2(z,80,fp); strcpy(cname[n],z);
  318.       fgets2(z,80,fp); strcpy(lname[n],z);
  319.       fgets2(z,80,fp); sens[n] = atosf(z);
  320.       fgets2(z,80,fp); Eauger[n] = atosf(z);
  321.       fgets2(z,80,fp); Ebind[n] = atosf(z);
  322.       fgets2(z,80,fp); Ratom[n] = atosf(z);
  323.    }
  324.    fclose(fp);
  325. }
  326.  
  327.  
  328. fgets2(s,n,fp)
  329. char *s;
  330. int n;
  331. FILE *fp;
  332. {
  333. int i,l;
  334.  
  335.    fgets(s,n,fp);
  336.    l = strlen(s);
  337.    for(i = 0; i <= l; i++) {
  338.       if(s[i] < 33) s[i] = 0;
  339.       if(s[i] == ';') s[i] = 0;
  340.    }
  341. }
  342.  
  343.  
  344.